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Abstract 

The magnetic perturbation due to the ferromagnetic test blanket modules (TBMs) may deteriorate fast ion confinement in ITER. 
This effect must be quantified by numerical studies in 3D. We have implemented a combined finite element method (FEM) - 
Biot-Savart law integrator method (BSLIM) to calculate the ITER 3D magnetic field and vector potential in detail. Unavoidable 
geometry simplifications changed the mass of the TBMs and ferritic inserts (FIs) up to 26%. This has been compensated for by 
modifying the nonlinear ferromagnetic material properties accordingly. Despite the simplifications, the computation geometry and 
the calculated fields are highly detailed. The combination of careful FEM mesh design and using BSLIM enables the use of the 
fields unsmoothed for particle orbit-following simulations. The magnetic field was found to agree with earlier calculations and 
revealed finer details. The vector potential is intended to serve as input for plasma shielding calculations. 
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1. Introduction 


The goal of the fusion reactor ITER is to demonstrate the tech¬ 
nological and scientific feasibility of fusion energy. The fol¬ 
lowing reactor, DEMO, has a mission to demonstrate the large- 
scale production of electrical power and tritium fuel self-suf¬ 
ficiency. One of the tasks of ITER is to be a testbed for DEMO 
components. ITER Test Blanket Modules (TBMs) will test the 
technology of tritium breeding modules for DEMO. Three of 
the eighteen ITER equatorial ports are reserved for these mod¬ 
ules. The material chosen for DEMO is ferritic steel (T). There¬ 
fore, the ITER TBMs will be made of ferromagnetic material 
that will get magnetized by the tokamak magnetic fields. The 
resulting local perturbation in the magnetic fields can deterio¬ 
rate the confinement of the plasma. Especially, the weakly col- 
lisional fast ions may find a local “hole” in the magnetic bottle 
and thus cause a hot spot on the first wall El. 

The TBM designs are checked for such threats by perform¬ 
ing simulations of fast ions 13j 01. The key ingredient for the 
simulations is a 3D magnetic field that includes the fields due 
to the magnetized components. This paper describes a method 


The used geometry describes the ferritic inserts and test blanket 
modules in detail We also calculate another useful quantity: 
the magnetic vector potential - an important input for related 
plasma response calculations. 

Input data. The key “raw” data we use to calculate the 3D 
field consists of the geometry and electric current in toroidal 
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field (TF) and poloidal field (PF) coils (including the central 
solenoid), the plasma equilibrium, and the geometry and mate¬ 
rial parameters of the ferromagnetic components. Two sets of 
components are considered: the TBMs and the ferritic inserts 
(FI) that are in place to reduce the variation of the toroidal field 
caused by having only 18 discrete TF coils. 

Methods. Our main tool is the commercial COMSOL Multi¬ 
physics finite element method (FEM) platform (ver. 4.4.0.195) 
and its AC/DC Module. For geometry simplifications we also 
used the SpaceClaim Engineer 3D direct modeler. In addition, 
MATLAB routines were used to prepare and deliver informa¬ 
tion to COMSOL as well as to build the COMSOL model. 

No FEM calculation can match the precision of a direct 
Biot-Savart law integration for magnetic fields from known coil 
geometry. To capitalise on this, we utilise a two-step COM¬ 
SOL calculation: First we perform a magnetization calculation. 
Then a follow-up permanent magnet COMSOL calculation ex¬ 
tracts the field due to the magnetization from the results of the 
magnetization calculation. The permanent magnet results from 
COMSOL are finally superimposed on the Biot-Savart law in¬ 
tegrated fields. These fields are produced with the recently ex¬ 
tended BioSaw code Q. 

In the magnetization calculation, the magnetic fields (mf) 
interface of the AC/DC module in COMSOL solves the vector 
potential A and the divergence control variable if/ from 

( V x H(|B|) = Je + V^ 

l V x A = B . (1) 

{ V-(A + V</r) = 0 

Symbol B denotes magnetic flux density and J e is the current 
density in coils or plasma. The nonlinear magnetic proper¬ 
ties are described by the function H(|B|). At the boundaries 
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of the computational domain, the boundary condition is set to 
n x A = 0, where n is the boundary’s normal vector. How¬ 
ever, the boundaries are effectively several kilometres from the 
tokamak, as will be explained in the context of equation [2] 

The permanent magnet model can be considered a reduced 
version of the above problem. We removed all coil and plasma 
currents and changed the constitutive relation inside the ferro¬ 
magnetic components to B = //o (H + M), where M is the mag¬ 
netization vector. COMSOL extracts the magnetization from 
the solution of the first step and effectively turns the ferromag¬ 
netic components into permanent magnets. The divergence con¬ 
dition variable is used for the permanent magnet model only 
when the vector potential needs to be evaluated. 

2. Geometry 

In this section we describe each major component present in the 
FEM model. The description justifies and explains all simplifi¬ 
cations made to the geometry before importing it to COMSOL. 

The 18 ITER toroidal field coils are D-shaped supercon¬ 
ducting coils enclosing the plasma. We reconstructed the 3D 
geometry by sweeping the coil cross section © over the oper¬ 
ating temperature spine curve Q. However, the curve was not 
perfectly smooth and continuous after it was drawn according 
to the specifications. Therefore a 2 nd degree Bezier curve (8| 
was fitted to a dense sampling of the CAD curve. The numer¬ 
ically smooth final curve is in a format natively supported by 
COMSOL. 

The poloidal field (PF) coils, the central solenoid (CS) and 
the plasma are assumed to be toroidally symmetric. We re¬ 
ceived the geometry information in EQDISK format, a quasi¬ 
standard in fusion community. ITER has 6 CS coils and 6 PF 
coils. The only change to this input was the elimination of a 5 
cm wide gap between the stacked CS coils by symmetrically 
stretching each coil vertically until they touched each other. 
Including this gap would have required small elements in the 
space between the coils without contributing significantly to the 
calculation. The current density was reduced to compensate for 
the increased cross-sectional area of the CS coils. The geome¬ 
try is shown in Fig. [T] 

The plasma cross-section was covered with a rectangular 
mesh (which was later swept toroidally) with the goal of ex¬ 
tending the mesh up to the plasma-facing components (Fig. 
[ljb)). The mesh nodes were carefully fitted to coincide with 
the grid used in the EQDISK calculation, though the pitch was 
a multiple of the EQDISK grid pitch. This grid setting aligned 
the element edges to cylindrical coordinates. This was useful, 
as the current densities were now always parallel to an element. 
Furthermore, the calculated fields would later be exported in the 
same coordinate system. 

The FI geometry we received represented the “configuration 
model” of the components, including the special elements at the 
NBI ports. This is a simplified description primarily intended 
for checking that no two components in the vast ITER CAD 
model collide. Therefore small details such as bolt holes had 
already been removed. Nevertheless, the CAD data was still 



Figure 1: (a) A single 20 degree sector including all the current-carrying do¬ 
mains. The calculation was performed in the complete 360 degree geometry, 
(b) The cross-sectional mesh of the toroidal current-carrying domains. The 
thick round curve is the plasma separatrix. 



Figure 2: Various 3D models for the ferromagnetic components, (a) a simple 
FI geometry for fast testing, (b) and (c) are the FI models used for calculations 
(no difference in results), (d) the original FI configuration model consisting of 
stacks of 40mm thick plates, which were merged in (c). The TBM geometry 
before (e) and after (f) simplification. 


far too detailed for finite element analysis (FEA) as a part of a 
model encompassing the whole ITER torus. 

We made several simplified FI models starting from the 
configuration model, that are shown in Fig. |2ja-d). The first 
simplification was to remove radial gaps. The FIs are made of 
stacks of 40 mm thick steel plates with 4 mm gaps between the 
plates. To expedite the calculations, several toroidal gaps less 
than 10 mm wide were removed and a few somewhat wider gaps 
were extended by 10 mm. Modifying the gaps did not change 
the final results. Many small surface details, such as bevels, 
were also removed to simplify the final model. 

Three equatorial ports of ITER are dedicated to TBMs, with 
room for two TBMs in each port. In our study, we placed 
a model of the European helium-cooled pebble bed (HCPB) 
TBM into all six TBM slots. In reality the different ITER part¬ 
ners will have their own test blanket module implementations 
in their slots. 

The CAD model for the HCPB is very complex, but we re¬ 
ceived a model which had cooling ducts and many other smaller 
features already removed. We further simplified the model by, 
e.g., removing piping, merging thin plates together and remov¬ 
ing air gaps. Figures |2je-f) show the geometry before and after 
simplification. 

The whole torus, including the coils, was enclosed inside a 
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so called “finite sphere” with a radius of 14 meters. This radius 
allowed us to fit the components inside the sphere with a mar¬ 
gin of several meters in most cases, but only about one meter 
near a PF coil. The same structure was also implemented in the 
permanent magnet model albeit with all the coils absent. The 
finite sphere was enclosed inside a spherical shell, dubbed the 
“infinite shell”. The name comes from the COMSOL feature 
we used to remap the radial coordinate g 2 = R 2 + z 2 inside the 
infinite shell in order to use boundary conditions at infinity: 


Removing small details changed the metal volumes of the 
TBMs and FIs. To compensate for this, we modified the mag¬ 
netic response of the materials, i.e., the H(|B|) function. We 
required that the simplifications would not change the magnetic 
moment M = f MdV of the objects at the uniform magnetic 
field limit. This resulted in the formula for a new H-B curve 
H*(B ), where the difference in metal volume is accounted for: 

H*(B) = H(B + {l- c}{B - g 0 H(B )}). (3) 


Q = Qo 


A Q 


Qo + A g — g 


( 2 ) 


where go is the inner radius of the infinite shell and A g is the 
thickness of the shell. In our model the outer surface of the shell 
was mapped to a distance of several kilometers. 


Here the volume ratio c is defined as c = Voriginai/ ^simplified, and 
it varies between 0.7355 and 0.738 for the different kinds of FI 
sectors and is 1.001 for the TBMs. 

4. Results 


3. Current densities and material properties for the Finite 
Element Analysis 

We included the following free currents in our model: the toroidal 
and poloidal field coils (including the central solenoid) and the 
toroidal plasma current. The magnitude of the current density 
was assumed to be uniform in each coil. For the circular PF 
and CS coils, the current direction is trivial to calculate, but for 
the D-shaped TF coil we assumed the direction to be parallel to 
the spine curve described in section [2] COMSOL has the func¬ 
tionality to create a curvilinear coordinate system within the TF 
coil, but as this seemed to cause numerical problems in our ge¬ 
ometry, we chose another route: a MATLAB routine returning 
a unit vector parallel to the spine curve of the TF coil was cre¬ 
ated. The magnitude of the toroidal plasma current density was 
calculated with a MATLAB routine from the current flux func¬ 
tion / and the pressure flux function p read from the EQDISK 
file using equation 3.3.6 in (9). (Note: there is a typo in that 
equation; po should be in the denominator.) 

All domains were assumed to have the same electrical prop¬ 
erties: relative permittivity s r = 1 and conductivity <x = 1 S/m. 
The current density was fixed in A/m 2 . Most domains had linear 
magnetic response with relative magnetic permeability p r = 1, 
while for the ferromagnetic components we set the magnetiza¬ 
tion M as a function of magnetizing field H by defining the H-B 
curve, or H(B). 

The FIs will be made of SS430 stainless steel, for which 
the B-H curve was the mean curve of the B-H curves in the ta¬ 
ble 4-1 of flOl . The magnetic properties data for the EUROFER 
steel ED of the TBMs is temperature dependent, but we made a 
conservative assumption of uniform 350°C. A two-piece linear 
model for the H-B curve was constructed from the three avail¬ 
able parameters. The first linear segment was assumed to pass 
through the origin, and the slope was calculated from the ra¬ 
tio of the coercive field H c and the remanent magnetization M r . 
The slope of the high field segment was vacuum permeability 
po, and the knee point was calculated by solving the location 
where the first segment passes through saturation magnetiza¬ 
tion M s . 


The COMSOL calculation produced 3D magnetic flux density 
B and vector potential A for various ITER scenarios. Figure 
[3] shows these fields for the 15 MA plasma current H-mode 
scenario during flat top phase. We then combined the COM¬ 
SOL results with Biot-Savart law integrated background fields 
and analysed the field. Figure [4] shows a toroidal field rip¬ 
ple map (a measure of FI performance), the homoclinic tangle 
near the X-point, and Poincare plots showing the structure in¬ 
side the plasma separatrix. Finally, the field was decomposed 
into toroidal Fourier components, which are shown in Figure [5] 
These form the input for spectral plasma response codes. The 
magnetic fields were verified against earlier work lfl2l and were 
found to agree quantitatively. Please check the supplementary 
material (available electronically) for illustration. 

5. Summary and Future work 

We have devised a method for calculating the magnetization of 
ITER ferromagnetic components using the finite element method, 
and combined the resulting magnetic flux density B and mag¬ 
netic vector potential A to Biot-Savart law integrated fields. The 
former can be used for studying, e.g., fast ion behaviour in the 
realistic ITER 3D field, while the latter provides the input for 
calculating the plasma response to the external perturbation. 

Future work. We are using the calculated fields in the ASCOT [13 ] 
suite of codes in order to simulate fast ion wall loads due to fu¬ 
sion alphas and heating neutral beams. 

The simulations are in progress, and we already reported 
that the detail level at which the FIs are modeled in the field cal¬ 
culation does not appear to play a significant role for fusion al¬ 
pha wall loads at least in the 15 MA H-mode case mm. The 
3D fields will also be delivered to our collaborators so that the 
plasma response to the external perturbations can be included in 
future wall load simulations. At the time of writing this article, 
we have not discovered strong changes in the wall loads due the 
European TBMs. Only a couple of ITER scenarios have so far 
been addressed. 
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Figure 3: The magnetic flux density B (T) and the total magnetic vector poten¬ 
tial A (T/m) as calculated by COMSOL. The total field is shown on the top row 
and the field due to the ferromagnetic components on the lower row. All im¬ 
ages show three orthogonal cut planes displaying three orthogonal components 
of the field, as indicated in the figure. The thin black lines indicate component 
and domain boundaries in the model. There are harmless numerical artifacts 
visible in A in the “infinite shell”, caused by remapping of the spatial coordi¬ 
nates. 
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Figure 5: Toroidal Fourier decomposition of the toroidal components of B and 
A fields. A single component on a poloidal plane is shown in figures (a), (b), 
(d) and (e). Figures (c) and (f) show the amplitude of the first 65 modes along 
the separatrix (black line). 
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Figure 4: (a) Poloidal and (b) toroidal Poincare plot showing the induced island 
structure within the plasma, (c) toroidal field ripple map, (d) the homoclinic 
tangle near X-point. 
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